home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C/C++ Users Group Library 1996 July
/
C-C++ Users Group Library July 1996.iso
/
listings
/
v_10_12
/
1012084a
< prev
next >
Wrap
Text File
|
1992-10-11
|
2KB
|
67 lines
/* bc program to calculate Chebyshef economized polynomial
** for evaluation of sin(x) */
/* use bc -l < listing1.bc to get s() function
** replace t() with c() to fit cos(x)
** a(x)/x for atan() */
define t(x) { /* sin(x)/x */
if(x==0)return(1.); /* l'Hopital's rule */
return (s(x) / x); /* put function to be fit here */ }
define b(x) { /* abs(x) */
if (x < 0) return (-x);
return (x); }
define m(x, y) { /* max(x,y) */
if (x > y) return (x);
return (y); }
n = 22; /* number of Chebyshef terms - need at least 3 more than
** are used below */
scale = 40;
/* set precision for dc -must be more than desired result */
p = a(1.) * 4; /* pi */
b = p * .25;
/* upper end of curve fit interval; b = 7/16 for atan() */
a = -b; /* lower end of interval */
/* chebft adapted from Press Flannery et al */
/* "Numerical Recipes" FORTRAN version */
for (k = 1; k <= n; ++k) {
c[k] = 0;
f[k] = t(c((k - .5) * p / n) * (b - a) * .5 + (b + a) * .5);
}
/* because of symmetry, even c[] are 0 */
for (j = 1; j <= n; j += 2) {
s = 0;
q = (j - 1) * p / n;
for (k = 1; k <= n; ++k) s += c(q * (k - .5)) * f[k];
/* old style bc requires =+ */
(c[j] = 2 / n * s); } /* display results */
/* skip even terms, which are 0 */
for (n = 5; n <= 19; n += 2) {
/* chebpc */
for (j = 1; j <= n; ++j) d[j] = e[j] = 0;
d[1] = c[n];
for (j = n - 1; j >= 2; --j) {
for (k = n - j + 1; 1 == 1; --k) {
s = e[k];
e[k] = d[k];
if(k==1)break;
d[k] = d[k - 1] * 2 - s; }
d[1] = c[j] - s; }
for (j = n; j >= 2; --j) d[j] = d[j - 1] - e[j];
d[1] = c[1] * .5 - e[1];
/* pcshft */
g = 2 / (b - a);
for (j = 2; j <= n; ++j) {
d[j] *= g; /* bc doesn't support multiple assignment */
g *= 2 / (b - a); }
for (j = 1; j < n; ++j) {
h = d[n];
for (k = n - 1; k >= j; --k) {
h = d[k] - (a + b) * .5 * h;
d[k] = h; }
}
"Chebyshev Sin fit |x|<Pi/4 coefficients"
" Maximum Rel Error:"
m(b(c[n + 2]), b(c[2])) / t(b);
for (i = 1; i <= n; i += 2) d[i];
}
WRAP_EOF